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Abstract 

In the paper, we propose a new calculation scheme for American options in the 
framework of a forward backward stochastic differential equation (FBSDE) . The well- 
known decomposition of an American option price with that of a European option 
of the same maturity and the remaining early exercise premium can be cast into the 
form of a decoupled non-linear FBSDE. We numerically solve the FBSDE by applying 
an interacting particle method recently proposed by Fujii & Takahashi (2012d), which 
allows one to perform a Monte Carlo simulation in a fully forward-looking manner. 
We perform the fourth-order analysis for the Black-Scholes (BS) model and the third- 
order analysis for the Heston model. The comparison to those obtained from existing 
tree algorithms shows the effectiveness of the particle method. 
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1 Introduction 



It has been known for a while that an American option value can be decomposed into 
that of the corresponding European option and an additional early exercise premium. 
Detailed discussions and other related references are available in Kim (1990) [23], Carr 
et.al. (1992) [7], Jacka (1991) [21], Rutkowski (1994) [32], Saito & Takahashi (2003) [33] 
as well as a textbook written by Karatzas & Shreve (1998) [23] . See also a recent work of 
Benth (2003) [2] which derives it from a dynamic programming approach. In this paper, we 
deal with a non-linear forward backward stochastic differential equation (FBSDE) obtained 
from the decomposition formula and calculate an American option price by solving it 
numerically. 

The framework of FBSDE was first introduced by Bismut (1973) [5], and then later 
extended by Pardoux and Peng (1990) [31] for general non-linear cases. Their financial ap- 
plications are discussed in details in, for example, El Karoui, Peng and Quenez (1997) [TT] 
and Ma and Yong (2000) [28]. There are increasing interests among researchers in FBS- 
DEs since their relevance for the analysis of various social phenomena is becoming more 
apparent in recent years. In fact, one can find FBSDEs in the valuation problem of the 
financial contracts in the presence of credit risk and/or funding cost of collaterals ( Duffie 
& Huang (1996) [10], Fujii & Takahashi (2012a) [13J, Crepey (2012) [8J, for examples. ). 
They are also relevant for the utility-indifference pricing in incomplete as well as con- 
strained markets ( Carmona (2009) [6j and references therein. ). In a recent book of 
Cvitanic and Zhang (2012) [9], the authors use FBSDEs to study the optimal contract 
theory in continuous time. 

Recently, Fujii & Takahashi (2012b) [2] has proposed a perturbative technique for 
generic non-linear FBSDEs. With the help of asymptotic expansion ( Takahashi (1999) [31] 
), it is possible to derive closed- form analytic expressions for both of the backward com- 
ponents. An explicit example for a quadratic-growth FBSDE appearing in the optimal 
portfolio problem in an incomplete market is available in Fujii & Takahashi (2012c) [15]. In 
the following paper, Fujii & Takahashi (2012d) [16] gave its numerical evaluation scheme 
based on an interacting particle method inspired by the work of McKean (1975) [29], 
which enables one to perform Monte Carlo simulation in a fully forward-looking man- 
ner. H The validity of its approximation is discussed recently by Takahashi & Yamada 
(2012b) [37J although it is still restricted to a decoupled non-linear setup. In the current 
paper, we apply this methodology to evaluate a non-linear FBSDE relevant for an Ameri- 
can option. Although there remains a small error when the option is far in the money, we 
shall see the effectiveness of the particle method in overall region. The current work not 
only gives a simple calculation scheme for American options but also serves as a concrete 
example showing the usefulness of the particle method to analyze non-linear FBSDEs and 
corresponding non-linear partial differential equations. 

It is closely related to the research with a long history on the branching Markov process and 
a certain class of semi-linear PDEs. For instance, see Fujita(1966) [T2], Ikeda, Nagasawa & Watan- 
abe(1965, 1966,1968) [17], [18], [19], Ikeda et.al.(1966, 1967) [20] and Nagasawa & Sirao (1969) [30]. 

2 A related but different approach was recently applied to evaluate CVA by Henry-Labordere (2012) |26] . 
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2 FBSDE formulation 



Let us take the probability space as (O, J 7 , Q), where Q is a risk- neutral probability mea- 
sure. We consider a generic process for the relevant stock price as 



dS t = (r t - yt)S t dt + S t a t ■ dW t , 



(2.1) 



where W is a d-dimensional Q-Brownian motion and J 7 is a natural filtration generated by 
W. All the stochastic processes are assumed to be J^-adapted. Here, r and y are processes 
for a risk-free interest rate and a dividend yield, respectively, a 6 W 1 is a volatility process. 

It is well-known (e.g. (2H [211 E21 [23]) that the price of an American option on S 
with a strike K and an expiry T can be expressed as 



f 3 u 1C ul{V u <V+(S u )} du 



where *f/ + (x) = max( x I'(x), 0) denotes a payoff function, which is 

Ct is a process denoting an instantaneous early exercise premium 

C t 

and 



ft 



(2.2) 



x — K for a Call 
K — x for a Put 



2/tiSt — rtK for a Call 
rtK — ytSt for a Put 



/3 t = exp( / r s ds 



(2.3) 



is a standard money-market account. 

In the remaining part of this section, we provide a simple heuristic derivation of 
Eq. f[2.2[) for completeness. Firstly, let us provide the decomposition principle of the 
Snell envelope for a continuous semimartingale. 



Proposition 1 Rutokowski (1994) 
Suppose X is a continuous semimartingale with canonical decomposition 

X = X + M + V 



(2.4) 



where Xq is a constant, M is a continuous local martingale with AIq = 0, and V denotes 
a continuous finite variation process with Vq = 0, whose decreasing component satisfies 
dV t d = v t dt for some adapted nonnegative process v. We assume that the condition 



E 



sup \X t \ 

0<t<T 



< CO 



satisfied. Let {T^}t^[o,T] be a family of {J- 1}- stopping times satisfying 



E[A>] = ess sup E[X r |J" t ], Vt € [0,T] . 

t<T<T 



(2.5) 



(2.6) 
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Then the following equation holds: 

E[X T *\T t ] = E[X T \T t ]-E 

Proof: See Appendix of [32J. 



JtT 



(2.7) 



For concreteness, let us choose a Call option with strike K as an example. We consider 
the dynamics of the discounted payoff process. By applying Ito formula, we obtain 



d{Pt\S t -K)- 



= -Pt l n{S t - K) + dt + ft" 1 \l {St > K} dSt + \s{X t - K)d(S) t \ 

= Pt ll {s t >K}St<y t ■ dW t 

+Pr 1 l {St >K}(rtK - y t S t )dt + \fr l 5{S t - K)d{S) u (2.8) 

where 5{-) is a Dirac delta function. More precisely speaking, the term involves the delta 
function is represented by the local time. For our intuitive derivation, however, the Dirac 
delta function is more useful to borrow a clear economic insight in a later stage. Now, 
applying Proposition Q] gives 



V t = ess sup /3 t E /3" 1 (S T -K)" 

t<r<T L 



Pt\St-K) 

J ^ l{r*=«} ^{Su^KyiVuSu ~ r u K)du - i^~ 1 *(5« - K)d(S) u ^ 



p-\S T -K)- 



+ m 



T 



1 {t * =u] P u 1 l {Su > K} (y u S u - r u K)di 



where, in the second equality, the last term vanishes due to the fact that the stock should 
be in-the-money region (S u > K) when the option is early exercised. It is now economically 
clear to see that the above result can be rewritten as 

cT 



v t = m PtHst-k) 



1 {Vu<{Su-K)+}Pu 1 {VuSu - r u K)du 



Note that 1 {V / u < (5ii _^) + }1{ 5ii> ^ } = l{v u <(s u -K)+} since the option value should always 
be positive. For more rigorous treatment, see the related proof in |32} [23] as well as [2]. 
The case for a Put option can be shown similarly. ■ 
Now, from Eq. (|2.2p . one can see 

Pt l Vt + I Pu l Cul{v u <*HSu)} du ( 2 - 9 ) 
Jo 

is a Q-martingale. Thus, we can conclude that the price of an American option satisfies 

dV t = r t V t dt - C t l {Vt < 9{St)} dt + Z t • dW t 

V T = ^ + {S T ) (2.10) 
dS t = (r t - yt)S t dt + S t a t ■ dW t , S = s 
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where Z E M rf is an appropriate J-f-adapted process that should be solved at the same 
time with V . It is a non-linear FBSDE with a decoupled dynamics of forward component, 
or the stock process S. Here, we have replaced \H r+ by \& in the indicator function since V 
should be clearly positive. In the next section, we carry out perturbative approximation 
procedures to solve the above FBSDE. 



3 Perturbative expansion and a particle method for FBSDE 

In [14j . a systematic approximation procedures for a generic non-linear FBSDE is given. 
It treats the non-linear driver of the FBSDE as a perturbation and converted the original 
system into a series of decoupled linear FBSDEs, for which the issue is equivalent to solve 
general European contingent claims. 

To apply the procedures, let us introduce perturbation parameter e as 



dV t (e) = r t V t (e) dt - eC t 6(*(S t ) - V t w )dt + Zf> ■ dW t 

v^ = *+(s T ) ( - Kl) 



where #(•) is the Heaviside step function. We now suppose that the solution of (|3.ip can 
be expanded as a power series of e: 

= V p + eVi W + e V/ 2) + e 3 v} 3) + • • • 

Economically speaking, we treat the early exercise premium as a perturbation and expand 
the price of American option around the corresponding European price. The method [2] 
allows to derive a series of linear FBSDEs specifying the dynamics of (V®, Z^')i>Q for each 
order of e. If the non-linear effects are sub-dominant and allow perturbative treatments, 
we can expect to obtain a reasonable approximation of the original model by setting e = 1 
at the end of the calculations. For the evaluation of an American option, the driver (or 
drift term) of the FBSDE is independent of the martingale component Z. Thus, in the 
following, we can focus on the level component V. 



3.1 Oth order 

In the Oth order, we have 



(0) _ m . V (°)A* , *(0) .„,, , 

(0) . T ._.„ . (3-2) 



dV t w = r t V t vu 'dt + Z^> ■ dW t 



which clearly represents the dynamics of the corresponding European option price. We 
can easily see that it is solved as 

V t (0) =p t E t [/3^ + (S T )\ Ft] ■ (3.3) 

Although there is no explicit expression of (|3.3p for a generic stock process, it is always 
possible to obtain its approximation by asymptotic expansion (See [34j [25j ESI [36] for the 
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details of asymptotic expansion.). It allows us, at least approximately, to have an explicit 
expression of as 

V t i0) = v^(t,X t ) (3.4) 

where X t = (St, r t , y t , a t , • • • ) contains all the relevant state processes. If necessary, ap- 
plication of Ito formula or using the process of Malliavin derivative (T>tXt) yields the 
corresponding martingale component Z^°\ 



3.2 1st order 

In the 1st order, the relevant FBSDE is given by 



dV, 



(1) 



r t V t (1) dt - C t e(V(S t ) - V W(t,X t ))dt + zf ] ■ dWt 



(1) 



,(1) 







(3.5) 



which is again linear and easy to integrate. We have 

rT 



v; 



(i) 



du/3 t E 



p- l c u e(^(s u )-v% 



(ay 



(3.6) 



where, denotes (u,X u ). is obtained by the similar arguments given in the 

previous subsection. Although it is possible to evaluate (|3.6p directly by Monte Carlo 
simulation, the time integration makes it rather time consuming. In fact, it soon becomes 
infeasible when one evaluates e-higher order expansion terms. 

In order to avoid the difficulty, we adopt an interacting particle method proposed in 
Fujii & Takahashi (2012d) [16]. We introduce an arbitrary J^-adapted strictly positive 
process {\t}t>o to define 



(i) 



exp 



(jf Xudu^Vp 



and 



Q,s = y ex P \ J ^uduj C s 
for s > t. Then, we have the SDE of vj\^ for the time component s (> t), 



(3.7) 
(3.8) 



dVK> = (r. + \ 8 )V£'ds - X s C t , s 9(^(S s 



't.s 

v {1) - 

v t,T — U 



,(°) 



)ds + ^du Z (i) . dWg 



(3.9) 



Since V$ = v/ 1 ^, we have 



V t 



(i) 



E 



J e-ft^+ x ^ du X s Ct, s 9(^(S s ) - vi 0) )ds 
l {T1< T}e- S * lr " du d t:T1 9(y(S. 



l {ri>t} 



(0)\ 



(3.10) 
(3.11) 



Here, t\ is a J^-stopping time associated with the first jump of Poisson process whose in- 
tensity process is given by {At}i>o- In contrast to (13.6f) . it is clear the expression of (13.111) 
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allows one-shot Monte Carlo simulation. More detailed explanation for Monte Carlo sim- 
ulation will be given in the later section. Although it is an interesting topic to obtain an 
optimal intensity process A that achieves the smallest variance in simulation, it is beyond 
the current scope of the paper. In the numerical examples, we simply use a constant 
intensity. 

Remark: In [16j . the intensity process A is assumed to be deterministic or an independent 
process for the other under lyings, which makes the evaluation of Malliavin derivatives re- 
quired for simpler. For the evaluation of American option, this assumption is not 
necessary since there is no need to obtain Z®. 



3.3 2nd order 

For the 2nd order case, the relevant equation is given by 



dV t 



(2) 



nV^dt + Ctd^iS, 



Vn 



(2) 



(0)V(1) 



V} 'dt + Z. 



(2) 



dW t 







(3.12) 



where 5(-) is a Dirac delta function as before. Since the FBSDE is linear, one can show 
easily that 

fT 



V 



(2) 



-A 



duE 



Ft 



(3.13) 



As mentioned in the previous section, the difficulty in a naive application of Monte Carlo 
simulation becomes much clearer now. At each point of time u £ [t, T] in a given path, 

one needs the value of Vu , which in turn requires to run Monte Carlo simulation as well 
as time integration. 

Therefore, let us define 



V, 



(2) 



t.s 



exp 



X u du)vP 



and use Ct, s as before. Then, for s > t, we have 

dvff = (r s + \ s )vffds + \ s d t>s 5($(S s ) - v(°A Vj^ds + e S' x « du Z^ ■ dW s 
with Vt t = 0. Thus, one obtains 



(3.14) 



(3.15) 



v; 



(2) 



y(2) 
v t,t 



-E 



"l{n>t} E 



l {ri<T}(- 



■JP r u du^ Ti§{Ti)v (l) 



Simple application of the tower property of iterated expectations gives 



V 



(2) 



Ft 



(3.16) 



(3.17) 



where t\ {t-i) is the first (second) jump time of the Poisson process with the intensity 



process A. Here, we have written S(t) = 5(9 (t) 
lighten the notations. 



and 8(t) = 0(tf(r) 



v^) to 
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3.4 3rd order 

In the 3rd order, the relevant dynamics becomes 



dV 



(3) 



r t V®dt + C t m*(S t )-v™)V t 



(0)W(2) 



= 0. 



-\d8(*{St) - vf ] ) {V t (l) ?}dt + Z t (3) • dW t , 



Here, the derivative of a Dirac delta function can be evaluated by approximating the delta 
function as a normal density function with a small variance, or using the integration-by- 
parts formula if possible. For the details of calculation, see the later sections treating 
numerical examples. After integration, we obtain 



V, 



(3) 



-Pt 



duE 



Ft 



T 



duE 



1 



Pu l C u -d5mS u )-VP){V^f 



Ft 



= exp 



Let us compress a convoluted expectation as before. Let us denote 

X u du\vj 3) 

and continue to use the simplified notations: 

(0), 



(3.18) 



(3.19) 



6(t) = 6mS t )-vn- 
Then, (|3.18|) is equivalent to 

dvff = (r s + X s )vffds + \ s C t , s {5(s)vW - ~d6(s)[V s W] 2 }ds + e^ x - du zf) ■ dW s 
with Vt t = 0, thus 



(3.20) 
(3.21) 



V t {3) = E 



l{n>t}K 



F t 
F t 



(3.22) 



Borrowing the idea from McKean [29] and use the tower property of iterated expectations, 
we finally obtain 



V 



+1 



(3) 



l{n>t} E 



-{ti<T2<T3< 



F 



{n>i} 



E 



p=i I 



(p) 



F t 
(3.23) 
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with the i-th jump time of the Poisson process denoted by Tj. 

In (|3.23p . p = {1,2} indicates one of the two particle groups. In both of the groups, 
the relevant state variables (or particles) follow the common diffusion dynamics ( those 
specified by BS or Heston models, for example) but those belong to different groups are 
independent each other, i.e., driven by the two independent set of Brownian motions. This 
particle representation compresses (E[ • I-Tvj) 2 into a single expectation operation. 

More concretely, for the evaluation of the second line, we use the branching diffusion 
method of McKean. For a each path of simulation, we 

(1) : update the diffusion process of the underlyings X = {S, r,y, a, • • • } in a standard way. 

(2) : do Poisson draw with intensity A at each step. 

(3) : if it draws a "jump" (or particles interact) at t\ < T, then the path yields the 
two identical copies of particles {X p } p= i 2 of the underlying states as its offspring, which 
continue to evolve according to the identical diffusion equations but driven by the two 
independent set of Brownian motions. 

(4) : for each particle group, we continue the Poisson draw of the second interaction until 
the maturity. 

(5) : finally, extract the following term: 

l^nle-^^dt^in) f[ l {Ti<T , <T} e-^ r ^C TuT r9(r^ (3.24) 
where is the second interaction time of each particle group. 

(6) : Repeat the procedures (1-5) and take the average of the values gathered in (5). 



3.5 4th order 

We can continue the expansion to an arbitrary higher order. In the 4th order, we have 



dV, 



M) 



i4 4) = o 



nV^dt + C t {k&smSt) - [V t {1) ] 3 - OWS) - [V®] 
+5(*(S t ) - « t (0) ) [V t {3) ] }dt + Z t (4) • dW t 



and hence 



(4) 



-Pt 
-Pt 
~Pt 



E 



1)13 



du 



T 



E 



du 



T 



E 



(3- l C u 5{nSu)-v^)V® 



du 



(3.25) 



Using the similar notations as in the previous sections, one can show that 
dV t f = (r s + X s )V t fds + X s C t ,s {^d 2 5(s)[VPf - d5(s)[VP][VP] + 5(s)V^\ ds 



+e f t s Kdu z (4) _ dWg 



(3.26) 
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and hence 



V t {A) = E 



^(r^du^ |_i. 5 2 <5(s)[y a) ] 3 + dS(s)[V s W}[VW} - 5(S)[VP}} ds 



Using the tower property and particle representation, the above result can be expanded 



as 



V, 



(i) 



p=i I 



(p) 



r T (p) 



x < 1 



>) 



{rKr^ <T 3 W <T} 



p=2 



T (P) <r| e-^ 3 ^ ds C iiT ( p) 5(T 2 (p) )C r ( p)iT ( p) 0(rf ) 



+ 1 {ri<T 2 <r 3 <r4<T}e" Jt 8 * Ct^^C^^ij^C^^^C^^Q^i) 

+ l{n<r 2 <T}^-^ 2 ^ rfs a i , Tl 5(r 1 )a ri , r2 ^(r 2 ) 



p=l 



2 



T (p) 
r 3 



ft 



4 Numerical Examples 

This section demonstrates the validity of our method proposed in the previous section 
through numerical experiments. 

4.1 Example 1: Black-Scholes model 

The first example is taken from Black-Scholes model: 



dSt/St = (r - y)dt + adW t , 



(4.1) 



where r, y and a are all nonnegative constants. We calculate the values up to the fourth 
order terms based on our scheme derived as f)3.11j) . ()3.17p . (|3.23p and (|3.2Tjl with 10 
million trials in Monte Carlo simulation. Here, we adopt the values reported in [22] as 
benchmarks. In particular, difficulty arises in differentiations up to the second order of 
the delta functions required for evaluation of the third (|3.23p as well as the fourth f|3.2Tj) 
order terms. Since the density function in Black-Scholes model is explicitly known, we 
are able to apply integration by parts (IBP) for computation of these terms in order to 
avoid differentiation of the delta functions. Moreover, we approximate each delta function 
by a normal density function with mean zero and a certain variance, which enables direct 
evaluation of the expectation. 



(3 
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Table 1: American Puts (T = 3,K = 100, a = 0.2, r = 0.08) 





So 


Benchmark 


0th 


1st 


2nd 


3rd 


4th 


y = 0.12 


80 


25.658 


24.777 


25.829 


25.854 


25.799 


25.739 




90 


20.083 


19.620 


20.174 


20.187 


20.158 


20.131 




100 


15.498 


15.252 


15.546 


15.553 


15.538 


15.518 




110 


11.803 


11.671 


11.830 


11.834 


11.826 


11.822 




120 


8.886 


8.814 


8.900 


8.902 


8.897 


8.894 


y = 0.08 


80 


22.205 


19.525 


23.553 


22.847 


22.265 


22.194 




90 


16.207 


14.676 


16.982 


16.522 


16.372 


16.316 




100 


11.704 


10.817 


12.151 


11.885 


11.800 


11.735 




110 


8.367 


7.847 


8.629 


8.473 


8.424 


8.409 




120 


5.930 


5.622 


6.081 


5.989 


5.962 


5.951 


y = 0.04 


80 


20.350 


14.589 


23.683 


22.236 


21.450 


20.348 




90 


13.497 


10.326 


16.120 


13.774 


13.573 


13.825 




100 


8.944 


7.168 


10.390 


9.132 


9.070 


8.706 




110 


5.912 


4.902 


6.720 


5.992 


5.957 


5.882 




120 


3.898 


3.315 


4.360 


3.953 


3.928 


3.827 


y = 0.00 


80 


20.000 


10.253 


24.338 


22.044 


20.892 


20.063 




90 


11.697 


6.783 


16.534 


12.950 


11.525 


11.959 




100 


6.932 


4.406 


9.590 


6.719 


7.004 


6.697 




110 


4.155 


2.826 


5.529 


4.066 


4.198 


4.286 




120 


2.510 


1.797 


3.232 


2.457 


2.506 


2.582 



Number of simulations = 10,000,000, A = 2, Number of time steps = 6000 



Error ratio 





So 


Benchmark 


0th 


1st 


2nd 


3rd 


4th 


y = 0.12 


80 


25.658 


-3.434% 


0.666% 


0.764% 


0.550% 


0.316% 




90 


20.083 


-2.305% 


0.453% 


0.518% 


0.373% 


0.239% 




100 


15.498 


-1.587% 


0.310% 


0.355% 


0.258% 


0.129% 




110 


11.803 


-1.118% 


0.229% 


0.263% 


0.195% 


0.161% 




120 


8.886 


-0.810% 


0.158% 


0.180% 


0.124% 


0.090% 


y = 0.08 


80 


22.205 


-12.069% 


6.071% 


2.891% 


0.270% 


-0.050% 




90 


16.207 


-9.447% 


4.782% 


1.944% 


1.018% 


0.673% 




100 


11.704 


-7.579% 


3.819% 


1.546% 


0.820% 


0.265% 




110 


8.367 


-6.215% 


3.131% 


1.267% 


0.681% 


0.502% 




120 


5.93 


-5.194% 


2.546% 


0.995% 


0.540% 


0.354% 


y = 0.04 


80 


20.35 


-28.310% 


16.378% 


9.268% 


5.405% 


-0.010% 




90 


13.497 


-23.494% 


19.434% 


2.052% 


0.563% 


2.430% 




100 


8.944 


-19.857% 


16.167% 


2.102% 


1.409% 


-2.661% 




110 


5.912 


-17.084% 


13.667% 


1.353% 


0.761% 


-0.507% 




120 


3.898 


-14.956% 


11.852% 


1.411% 


0.770% 


-1.821% 


y = 0.00 


80 


20 


-48.735% 


21.690% 


10.220% 


4.460% 


0.315% 




90 


11.697 


-42.011% 


41.352% 


10.712% 


-1.470% 


2.240% 




100 


6.932 


-36.440% 


38.344% 


-3.073% 


1.039% 


-3.390% 




110 


4.155 


-31.986% 


33.069% 


-2.142% 


1.035% 


3.153% 




120 


2.51 


-28.406% 


28.765% 


-2.112% 


-0.159% 


2.869% 



error ratio = 100*(value-benchmark)/benchmark 
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Table 2: American Calls (T = 3, K = 100) 





So 


Benchmark 


0th 


1st 


2nd 


3rd 


4th 


a 


= 


0.2 


80 


2.580 


2.241 


2.847 


2.612 


2.602 


2.566 


r 


= 


0.03 


90 


5.167 


4.355 


5.822 


5.240 


5.199 


5.140 


V 


= 


0.07 


100 


9.066 


7.386 


10.453 


9.204 


9.121 


8.877 








110 


14.443 


11.331 


17.036 


14.763 


14.566 


14.322 








120 


21.414 


16.117 


25.307 


23.173 


22.023 


21.044 


a 


= 


0.4 


80 


11.326 


10.309 


11.998 


11.475 


11.399 


11.373 


r 


= 


0.03 


90 


15.722 


14.162 


16.769 


15.975 


15.858 


15.770 


y 




0.07 


100 


20.793 


18.532 


22.318 


21.132 


20.948 


20.851 








110 


26.494 


23.363 


28.609 


26.902 


26.672 


26.536 








120 


32.781 


28.598 


35.599 


33.390 


33.000 


32.681 


a 




0.3 


80 


5.518 


4.644 


6.254 


5.564 


5.562 


5.511 


r 




0.00 


90 


8.842 


7.269 


10.197 


8.951 


8.914 


8.763 


V 




0.07 


100 


13.142 


10.542 


15.407 


13.285 


13.178 


12.989 








110 


18.453 


14.430 


22.004 


18.655 


18.726 


18.144 








120 


24.791 


18.882 


30.019 


25.587 


23.554 


24.215 


a 




0.3 


80 


12.146 


12.133 


12.148 


12.148 


12.147 


12.146 


r 




0.07 


90 


17.368 


17.343 


17.373 


17.374 


17.372 


17.372 


V 




0.03 


100 


23.348 


23.301 


23.359 


23.360 


23.357 


23.355 








110 


29.964 


29.882 


29.980 


29.982 


29.977 


29.976 








120 


37.104 


36.972 


37.130 


37.134 


37.125 


37.120 



Number of simulations = 10,000,000, A = 2, Number of time steps = 6000 



Error ratio 





So 


Benchmark 


0th 


1st 


2nd 


3rd 


4th 


a 


= 0.2 


80 


2.58 


-13.140% 


10.349% 


1.240% 


0.853% 


-0.543% 


r 


= 0.03 


90 


5.167 


-15.715% 


12.677% 


1.413% 


0.619% 


-0.523% 


V 


= 0.07 


100 


9.066 


-18.531% 


15.299% 


1.522% 


0.607% 


-2.085% 






110 


14.443 


-21.547% 


17.953% 


2.216% 


0.852% 


-0.838% 






120 


21.414 


-24.736% 


18.180% 


8.214% 


2.844% 


-1.728% 


a 


= 0.4 


80 


11.326 


-8.979% 


5.933% 


1.316% 


0.645% 


0.415% 


r 


= 0.03 


90 


15.722 


-9.922% 


6.659% 


1.609% 


0.865% 


0.305% 


V 


= 0.07 


100 


20.793 


-10.874% 


7.334% 


1.630% 


0.745% 


0.279% 






110 


26.494 


-11.818% 


7.983% 


1.540% 


0.672% 


0.159% 






120 


32.781 


-12.760% 


8.596% 


1.858% 


0.668% 


-0.305% 


a 


= 0.3 


80 


5.518 


-15.839% 


13.338% 


0.834% 


0.797% 


-0.127% 


r 


= 0.00 


90 


8.842 


-17.790% 


15.325% 


1.233% 


0.814% 


-0.893% 


y 


= 0.07 


100 


13.142 


-19.784% 


17.235% 


1.088% 


0.274% 


-1.164% 






110 


18.453 


-21.801% 


19.243% 


1.095% 


1.479% 


-1.675% 






120 


24.791 


-23.835% 


21.088% 


3.211% 


-4.990% 


-2.323% 


a 


= 0.3 


80 


12.146 


-0.107% 


0.016% 


0.016% 


0.008% 


0.000% 


r 


= 0.07 


90 


17.368 


-0.144% 


0.029% 


0.035% 


0.023% 


0.023% 


V 


= 0.03 


100 


23.348 


-0.201% 


0.047% 


0.051% 


0.039% 


0.030% 






110 


29.964 


-0.274% 


0.053% 


0.060% 


0.043% 


0.040% 






120 


37.104 


-0.356% 


0.070% 


0.081% 


0.057% 


0.043% 



error ratio = 100*(value-benchmark)/benchmark 
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±auie o 


: American Calls (T = 0.5 


, & — 


iuu, r — 


n m ii 

U-Uc>, y — 


0.0 1 ) 




So 


Benchmark 


0th 


ist 




ora 


4tn 


cr = 0.2 


80 


0.219 


0.215 


0.222 


0.220 


0.220 


0.220 




90 


1.386 


1.345 


1.413 


1.391 


1.389 


1.389 




100 


4.783 


4.578 


4.920 


4.807 


4.795 


4.791 




110 


11.098 


10.421 


11.569 


11.172 


11.137 


11.221 




120 


20.000 


18.302 


20.536 


20.350 


20.238 


20.092 


cr = 0.4 


80 


2.689 


2.651 


2.710 


2.695 


2.692 


2.693 




90 


5.722 


5.622 


5.778 


5.736 


5.729 


5.730 




100 


10.239 


10.021 


10.365 


10.272 


10.257 


10.261 




110 


16.181 


15.768 


16.424 


16.241 


16.214 


16.232 




120 


23.360 


22.650 


23.778 


23.459 


23.411 


23.395 



Number of simulations = 10,000,000, A = 2, Number of time steps = 1000 



Error ratio 





So 


Benchmark 


0th 


1st 


2nd 


3rd 


4th 


a = 0.2 


80 


0.219 




-1.826% 


1.370% 


0.457% 


0.457% 


0.457% 




90 


1.386 




-2.958% 


1.948% 


0.361% 


0.216% 


0.216% 




100 


4.783 




-4.286% 


2.864% 


0.502% 


0.251% 


0.167% 




110 


11.098 




-6.100% 


4.244% 


0.667% 


0.351% 


1.108% 




120 


20 




-8.490% 


2.680% 


1.750% 


1.190% 


0.460% 


a = 0.4 


80 


2.689 




-1.413% 


0.781% 


0.223% 


0.112% 


0.149% 




90 


5.722 




-1.748% 


0.979% 


0.245% 


0.122% 


0.140% 




100 


10.239 




-2.129% 


1.231% 


0.322% 


0.176% 


0.215% 




110 


16.181 




-2.552% 


1.502% 


0.371% 


0.204% 


0.315% 




120 


23.36 




-3.039% 


1.789% 


0.424% 


0.218% 


0.150% 



error ratio = 100* (value-benchmark) /benchmark 



As for the third order term, consulting the results based on IBP, we are capable of 
determining an appropriate size of the variance for each normal density applied in the 
approximation of a delta function. Unfortunately, however, this IBP method does not 
yield stable results for some cases in computing the fourth-order term. It is clear that 
we want to use a small enough variance for the normal density so that it is a reasonable 
approximation of the delta function. On the other hand, too small variance increases the 
variation (dispersion) of simulation result. Therefore, we change the variance from some 
large value to a smaller one gradually for a given number of simulation paths and picks 
up the smallest value beyond which the variation (dispersion) starts to increase. 

This scheme can be applied to general cases where the density functions of the under- 
lying models are not explicitly available. In fact, we adopt this approach for the numerical 
example for the Heston model in the next subsection H. Of course, there is no guarantee 
that the choice of variance that gives the smallest dispersion in simulation also yields 
the smallest bias in the numerical result. However, the numerical result suggests that 
the method produces accurate enough approximation for practical use given a reasonable 

3 Notice that there is no need to use unnecessarily small variance for the approximation of delta function. 
Intuitively speaking, the delta function within the expectation operation extracts the density where its 
argument vanishes. Thus, as long as the density functions of the underlyings do not change significantly 
within a given range, one can use it as a variance of the normal density function as an approximation of 
the corresponding delta function. 
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number of simulation paths. Note here that the large number of paths used in this example 
is to confirm the convergence of higher order expansions. For practical pricing purpose, 
as can be seen in the following example of Heston model, there is no need to run such a 
large number of simulation trials. 

Table [T] presents the result for American put options with T = 3, K = 100, a = 0.2 and 
r = 0.08, which confirms that the error ratios become improved in the results up to the 
third or the fourth order comparing with those up to the first and the second orders. In 
total, the approximations up to the fourth order provide the most precise ones in terms of 
the error ratios. Note also that for the dividend rate y = 0.12 and 0.07, adding the fourth 
order term to the third one makes the accuracies of the approximations improved, while 
for y = 0.04 and 0.00, it makes the accuracies worse in three and four out of the five cases, 
respectively. Table [2] (T = 3, K = 100) and Table H (T = 0.5, K = 100, r = 0.03, y = 0.07) 
present the results for American call options, which shows the error ratios become smaller 
by adding the third or /and fourth order terms. 

4.2 Example 2: Heston model 

The next example takes Heston model (I4.2p : 

dS t = (r-y)Stdt + StVFfdW u (4.2) 
dv t = £(0 - u t )dt + T}y/u t (j>dWx,t + y/l- p 2 dW 2 ,t) ; ^ = a 2 , (4.3) 

where £, r\ and 6 are positive constants such that ^rj > 9 2 /2 and WitJLWsj,*- Then, 
we compute the approximate values for American put prices (T = 0.25, 0.5, K = 100, 
r = 0.05, y = 0, rj = 0.1, £ = 3.0, 9 = 0.04) up to the third order based on our scheme, 
(I3~TT]1 . (|3TT71) and (^23]) with 50,000 trials in Monte Carlo simulation. Here, we adopt [I] 
as the benchmark values, in which a two-dimensional tree with two hundred time steps 
and a Control Variate technique is applied. Moreover, an asymptotic expansion method, 
particularly, the equation appearing in p. 113 of [36 is used for computation of European 
option prices, that is, 0-th order ii' ', 

Table H] demonstrates that our method works effectively in the Heston model, which 
suggests its applicability to the pricing problem of American options for other multi- 
dimensional models, too. The numerical result shows that the expansion up to the third 
order improves the accuracies in most of the cases. For the choice of the normal density 
as a approximate delta function, we have used the variance found to work well in the 
previous BS model. In the case where it produces too much dispersion in simulation, we 
have applied the general methodology explained in the previous subsection to pick up an 
appropriate variance. The example shows that the relatively small number of simulation 
trials is enough to obtain a reasonable accuracy for the practical use. In Table [5j we have 
given the numerical results with larger number of simulations 500, 000 for the same set of 
American options with T = 0.5 in Table HI Although the improvement of accuracy from 
the second to the third order approximation becomes more robust in this case, one can 
observe that the size of the change in option prices is rather small. 

Remark: 

Although higher order integration is required, the direct evaluation of (|3. 13|) ( and corre- 
sponding expressions in other orders ) is also possible once we know the transition density 
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of the underlying states. For the diffusion models, it is always possible to obtain approx- 
imation using the asymptotic expansion [34J. If there exists efficient enough integration 
technique, such as Gaussian quadrature and its extension, it could provide another pricing 
technique. In fact, in BS model, we have compared this semi-analytic results (by brute 
force integration within ±5-sigma range and using a normal density function with variance 
of lbp of the stock process at each time as an approximation for the delta function) to 
those obtained from the particle method up to the second order terms. We confirmed the 
consistency between their numerical results. 

Table 4: American Puts in Heston model (K = 100, r = 0.05, y = 0,ry = 0.1, £ = 3.0,6 = 
0.04) 



T = 0.25 





So 


Benchmark 


0th 


1st 


2nd 


3rd 


ER(Oth) 


ER(lst) 


ER(2nd) 


ER(3rd) 


p = —0.1 

cr = 0.2 


90 
100 
110 


10.171 
3.475 
0.774 


9.643 
3.374 
0.758 


10.507 
3.556 
0.783 


10.293 
3.486 
0.775 


10.141 

3.481 
0.775 


-5.18 % 
-2.89 % 
-2.03 % 


3.31 % 
2.35 % 
1.21 % 


1.20 % 
0.34 % 
0.18 % 


-0.29 % 
0.19 % 
0.14 % 


p = -0.7 

cr = 0.2 


90 
100 
110 


10.121 
3.481 
0.842 


9.573 
3.383 
0.829 


10.455 
3.559 
0.854 


10.253 
3.493 
0.845 


10.101 

3.487 
0.845 


-5.41 % 
-2.81 % 
-1.53 % 


3.31 % 
2.24 % 
1.43 % 


1.31 % 
0.36 % 
0.37 % 


-0.19 % 
0.18 % 
0.41 % 


p = -0.1 

cr = 0.4 


90 
100 
110 


12.182 
6.496 
3.091 


11.896 
6.379 
3.047 


12.347 
6.572 
3.118 


12.190 
6.504 
3.092 


12.173 
6.501 
3.092 


-2.35 % 
-1.80 % 
-1.43 % 


1.36 % 
1.17 % 
0.85 % 


0.07 % 
0.12 % 
0.03 % 


-0.08 % 
0.08 % 
0.02 % 


p = -0.7 

cr = 0.4 


90 
100 
110 


12.112 
6.490 
3.146 


11.832 
6.377 
3.104 


12.291 
6.565 
3.180 


12.132 
6.505 
3.157 


12.116 
6.503 
3.157 


-2.31 % 
-1.74 % 
-1.31 % 


1.47 % 
1.16 % 
1.11 % 


0.16 % 
0.23 % 
0.36 % 


0.03 % 
0.20 % 
0.37 % 






Number of simulations = 50,000, A = 4, Number of time steps 
ER = 100* (value-benchmark) /benchmark 

T = 0.5 


= 2000 








So 


Benchmark 


0th 


1st 


2nd 


3rd 


ER(Oth) 


ER(lst) 


ER(2nd) 


ER(3rd) 


p = -0.1 

cr = 0.2 


90 
100 
110 


10.648 
4.647 
1.683 


9.864 
4.423 
1.624 


11.236 
4.835 
1.733 


10.752 
4.676 
1.692 


10.532 
4.665 
1.693 


-7.36 % 
-4.83 % 
-3.50 % 


5.53 % 
4.04 % 
2.94 % 


0.98 % 
0.61 % 
0.55 % 


1.09 % 
0.38 % 
0.55 % 


p = -0.7 

cr = 0.2 


90 
100 
110 


10.564 
4.664 
1.787 


9.766 
4.443 
1.732 


11.183 
4.844 
1.837 


10.688 
4.684 
1.798 


10.490 
4.678 
1.797 


-7.55 % 
-4.73 % 
-3.08 % 


5.87 % 

3.88 % 
2.79 % 


1.18 % 
0.43 % 
0.61 % 


0.70 % 
0.31 % 
0.52 % 


p = -0.1 

cr = 0.4 


90 
100 
110 


13.314 
8.008 
4.545 


12.712 
7.705 
4.399 


13.664 
8.207 
4.642 


13.375 
8.070 
4.567 


13.283 
8.021 
4.550 


-4.52 % 
-3.78 % 
-3.21 % 


2.63 % 
2.48 % 
2.12 % 


0.46 % 
0.77 % 
0.48 % 


0.23 % 
0.16 % 
0.09 % 


p = -0.7 

cr = 0.4 


90 
100 
110 


13.217 
8.000 
4.620 


12.625 
7.705 
4.479 


13.602 
8.196 
4.709 


13.314 
8.048 
4.627 


13.229 
8.012 
4.612 


-4.48 % 
-3.69 % 
-3.04 % 


2.91 % 
2.46 % 
1.93 % 


0.73 % 
0.60 % 
0.16 % 


0.09 % 
0.15 % 
0.17 % 



Number of simulations = 50,000, A = 8, Number of time steps = 2000 
ER = 100* (value-benchmark) /benchmark 



5 Conclusions 

This paper proposed a new calculation technique for American options in an FBSDE 
framework. The well-known decomposition of an American option price with that of the 
corresponding European option and additional early exercise premium can be written in a 
form of a decoupled non-linear FBSDE. We have used the recently proposed perturbation 
technique of FBSDE with an interacting particle method to obtain numerical results. We 
have tested the effectiveness of our approximation by comparing the numerical results 
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Table 5: The same setup with T = 0.5 in Table H] but using larger number of simulation. 



r = 0.5 





So 


Benchmark 


0th 


1st 


2nd 


3rd 


ER(Oth) 


ER(lst) 


ER(2nd) 


ER(3rd) 


p = 


-0.1 


90 


10.648 


9.864 


11.259 


10.758 


10.540 


-7.36 % 


5.74 % 


1.03 % 


-1.01 % 


a = 


0.2 


100 


4.647 


4.423 


4.831 


4.674 


4.653 


-4.83 % 


3.95 % 


0.57 % 


0.11 % 






110 


1.683 


1.624 


1.729 


1.688 


1.683 


-3.50 % 


2.70 % 


0.28 % 


-0.02 % 


P = 


-0.7 


90 


10.564 


9.766 


11.173 


10.696 


10.457 


-7.55 % 


5.77 % 


1.26 % 


-1.01 % 


a = 


0.2 


100 


4.664 


4.443 


4.841 


4.688 


4.672 


-4.73 % 


3.81 % 


0.53 % 


0.19 % 






110 


1.787 


1.732 


1.839 


1.800 


1.795 


-3.08 % 


2.86 % 


0.71 % 


0.44 % 


P = 


-0.1 


90 


13.314 


12.712 


13.676 


13.384 


13.311 


-4.52 % 


2.72 % 


0.53 % 


-0.02 % 


a = 


0.4 


100 


8.008 


7.705 


8.202 


8.043 


8.007 


-3.78 % 


2.41 % 


0.44 % 


-0.01 % 






110 


4.545 


4.399 


4.643 


4.561 


4.546 


-3.21 % 


2.14 % 


0.34 % 


0.02 % 


P = 


-0.7 


90 


13.217 


12.625 


13.582 


13.292 


13.216 


-4.48 % 


2.76 % 


0.57 % 


-0.01 % 


a = 


0.4 


100 


8.000 


7.705 


8.194 


8.039 


8.003 


-3.69 % 


2.42 % 


0.49 % 


0.04 % 






110 


4.620 


4.479 


4.718 


4.640 


4.625 


-3.04 % 


2.12 % 


0.43 % 


0.10 % 



Number of simulations = 500,000, A = 8, Number of time steps = 2000 
ER = 100* (value-benchmark) /benchmark 



to those obtained from existing tree algorithms. Although there remains some subtlety 
for choosing an appropriate variance for the normal density function as a proxy of the 
Dirac delta function, the proposed method for the variance choice yields accurate enough 
approximations for BS as well as Heston models. In the paper, we could only test a 
narrow range of parameters with relatively short expiries of options due to the limitation 
of existing benchmark results. However, the results are quite encouraging to suggest that 
our perturbation technique combined with an interacting particle method can be applied 
to much broader range of models and parameters. 
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